;; Utilities for the non-robust linear algebra functions
;; Copyright (C) 2025  Rocks Mazama
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <https://www.gnu.org/licenses/>.

(ns mazama.linalg.utilities
  (:require
   [mazama.linalg.matrix :as matrix]
   [mazama.linalg.sole :as sole]))

(defn +v
  "Add two vectors, `u` and `v`.
  If the number of components differ then the result is undefined."
  ([u] u)
  ([u v] (map + u v)))

(defn *v
  "Multiply a scalar `c` and a vector `v`."
  [c v]
  (map #(* % c) v))

(defn linear-combination
  "Create a linear combination from the vectors `vs` and coefficients `cs`.
  If the number of elements or components differ then the result is undefined."
  [vs cs]
  (let [product (map *v cs vs)]
    (if (empty? product)
      []
      (reduce +v product))))

(defn dot-product
  "Compute the dot product of two vectors.
  If the number of components differ then the result is undefined."
  [u v]
  ;; Slower alternative...
  ;; (reduce + 0 (map * u v))
  (loop [u u
         v v
         result 0]
    (if (or (empty? u) (empty? v))
      result
      (recur (rest u) (rest v) (+ result (* (first u) (first v)))))))

(defn *mv
  "Multiply a matrix A and a vector `x` to form the vector Ax."
  [{:keys [ncols data]} x]
  (->> data
       (map * (cycle x))
       (partition ncols)
       (map #(reduce + %))))

(defn *mm
  "Multiply a matrix `A` by a matrix `B` to form the matrix AB."
  [A B]
  ;; Slower alternative...
  ;;  (->> B
  ;;       (matrix/columns)
  ;;       (map #(*mv A %))
  ;;       (matrix/columns->matrix))
  (let [bcolumns (matrix/columns B)]
    (loop [arows (matrix/rows A)
           bcols bcolumns
           data []]
      (cond
        (empty? arows) {:ncols (:ncols B) :data data}
        (empty? bcols) (recur (rest arows) bcolumns data)
        :else (recur arows
                     (rest bcols)
                     (conj data (dot-product (first arows) (first bcols))))))))

(defn linear-combination?
  "Is `v` a linear combination of the vectors `vs`?"
  [v vs]
  (sole/solve (matrix/columns->matrix vs) {:ncols 1 :data v}))

(defn linearly-independent?
  "Are the vectors in the collection `vs` linearly independent?"
  [vs]
  (let [A (matrix/columns->matrix vs)
        b {:ncols 1 :data (repeat (/ (count (:data A)) (:ncols A)) 0)}
        x (sole/solve A b)]
    (and x (nil? (:nullspace x)))))

(defn random-seq [prng count]
  (repeatedly count prng))

(defn prng-bound [seed bound]
  (let [frand (java.util.Random. seed)]
    (fn [] (.nextInt frand bound))))

(defn prng-min-max
  "Return a pseudo-random number generator using seed.
  Each call to the prng returns min <= value <= max."
  [seed min max]
  (let [bound (inc (- max min))
        rndm (prng-bound seed bound)]
    (fn [] (+ min (rndm)))))

(defn random-matrix
  "Return a nxm matrix using prng to generate values"
  [prng n m]
  {:ncols m :data (random-seq prng (* n m))})

(defn random-sole-unique-solution
  "Return a System of Linear Equations Ax=b that has a unique solution.
    A is an nxn matrix,
    x is an nxp matrix,
    b is an nxp matrix.
    The elements of A and x are generated by prng."
  [prng n p]
  (let [A (random-matrix prng n n)
        x (random-matrix prng n p)
        b (*mm A x)]
    (vector A x b)))
